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reference  fluid  static  pressure 
body  radius 

geometric  parameter,  defined  in  Eq.  (22) 
velocity  components  in  the  x,  y,  z-directions 
transformed  variables,  defined  in  Eq.  (25) 
total  velocity 

shear  velocity  (=  V~\Jq  I 
components  of  Reynolds  stress 
nonorthorgonal  coordinates 
Cartesian  coordinates 
x  value  at  separation 
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angle  between  x  and  /  coord i nates 
geometric  parameters,  defined  in  Eq.  (6) 
an  intermittency  factor  in  transition  region,  defined  in  (llj 
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ABSTRACT 


Thick  boundary-layer  equations,  including  the  variation  of  the  metric 
coefficients  with  distance  normal  to  the  surface,  are  solved  from  the  nose 
to  the  tail  of  an  elongated  three  dimensional  body.  Iteration  between 
potential  and  boundary-layer  flow  solutions  is  adopted  to  obtain  a 
converged  solution  in  the  stern  region.  The  complete  procedure  makes  use 
of  Keller’s  box  numerical  scheme  in  a  non-orthogonal  curvilinear 
coordinate  system  and  an  isotropic  eddy -viscosity  model  for  the  Reynolds 
stress  tensor  in  both  streamwise  and  crosswise  directions.  The  numerical 
method  allows  the  calculation  of  flows  in  which  the  crosswise  component 
of  velocity  contains  regions  of  flow  reversal  across  the  boundary  layer.  The 
in  viscid  pressure  is  determined  by  using  the  Hess  Smith  method.  The 
program  is  validated  for  an  axisymmetric  bodv  and  then  the  calculations 
are  compared  with  measurements  for  two  three  dimensional  bodies  having 
elliptical  transverse  sections. 
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INTRODUCTION 

When  compared  to  the  thin  hull  boundary  laser,  much  less  progress  has  been  made  in  the 
prediction  of  the  three-dimensional  (3  D)  thick  stern  turbulent  boundars  laser  and  wake  which  invoke 
higher-order  effects.  Only  three  out  of  the  seventeen  prediction  methods  in  the  SSPA  ITTC  1980  Ship 
Boundary  Layer  Workshop1 ’considered  higher order  effects  This  is  because  of  (ai  the  complexity  of  the 
3  D  llow.  and  (bl  the  existence  of  flow  separation  and  trailing  vorticits  Complete  computational  solutions 
of  the  full  Navier  Stokes  equations  are  still  limited  to  laminar  time  averaged  turbulent  Hows.  These 
computations  require  the  use  of  computers  basing  large  cure  memories  and  high  speed  processors.  Seeking 
solutions  for  higher  order  equations,  whose  complexits  lies  somewhere  between  that  ot  the  Resnolds  and 
the  thin  boundars  layer  equations  for  turbulent  flow,  is  a  current  ongoing  research  topic  in 
computational  ship  hydrodynamics. 

Two  different  approaches  exist  for  solving  higher  order  equations  I  he  brsi  approach  simplifies  the 
Reynolds  equations  by  discarding  the  terms  related  to  streamwise  diffusion  1  he  resulting  equations  are 
called  the  partially  parabolic  equations.  1. samples  of  this  approach  are  contained  in  References  2.  3.  and 
4.  The  second  approach  generalizes  the  first  order  boundars  laser  equations  bv  introducing  terms  related 
to  the  second order  curvature  effects.  The  complete  solution  procedure  includes  viscous  mviscid  inter 
action  to  account  for  deficiencies  in  the  governing  equations  as  well  as  m  the  numerical  solution 
techniques.  This  approach  has  been  widcls  used  bv  the  aircraft  mdustrv.  applications  to  ship  flows  are 

*A  complete  listing  of  references  is  given  on  page  35. 


given  in  References  5  and  6.  The  parabolic  nature  of  the  equations  comprising  the  latter  approach 
indicates  the  downstream  pressure  field  does  not  influence  upstream  flow  characteristics.  Through  the 
viscid  inviscid  interaction  process,  the  solutions  are  corrected  by  successive  global  iterations.  Although 
Abdelmequid  et  al.2  has  made  some  comparisons  between  the  two  approaches,  he  unfortunately  adopted 
the  second  approach  without  including  interactions.  It  is  still  not  known  which  approach  gives  better 
predictions. 

In  this  paper,  the  numerical  method  of  Cebeci  et  al.7  for  solving  3-D  boundary  layers  is  expanded 
to  include  the  curvature  terms  for  thick  stern  boundary  layers.  Several  new  geometrical  parameters  are 
introduced  by  using  the  original  nonorthogonai  coordinate  system.  The  boundary  conditions  at  the  edge 
of  the  thick  boundary  layer  are  updated  in  the  global  iterations  between  viscous  and  inviscid  solutions. 
Some  encouraging  results  have  been  obtained  for  an  axisymmetric  flow  and  two  3-D  flows  with 
separation.  In  their  survey  of  boundary-layer  theory  for  flow  with  separation.  McDonald  and  Briley8 
concluded  that  2  D  flows  with  mild  separation  can  be  treated  by  applying  an  interaction  procedure  to  the 
boundary  layer  equations.  The  present  results  shed  some  light  on  extending  this  approach  to  3-D  flows. 

COORDINATE  SYSTEM 

A  curvilinear  nonorthogonai  surface  fitted  coordinate  system  is  used  to  represent  ship  hull  surfaces. 
Figure  I  depicts  this  coordinate  system  (x,  y,  /.)  along  with  Cartesian  coordinates  lx,  y.  z).  Here,  a  is  the 
angle  between  x  and  z,  and  y  is  normal  to  the  hull  surface.  The  potential  flow  solution  for  a  ship  hull  is 
solved  in  terms  of  the  Cartesian  coordinates  and  the  boundary-layer  equations  are  expressed  in  the  (x.  y, 
/)  coordinate  system.  We  select  x  =  x,  and  z  lies  in  the  yz-plane.  Because  the  coordinate  system  is 
nonorthogonai.  we  are  free  to  select  the  value  of  z  in  the  plane  to  satisfy  the  condition  that  the 
boundary  line  of  the  ship  hull  is  coincident  with  the  /  =  constant  coordinate  line.  In  this  paper,  z  is 
determined  by  mapping  each  yz  crossplane  into  a  half  unit  circle  by  the  numerical  method  developed  by 
Halsey1'.  The  polar  angle,  normalized  by  n  on  the  unit  circle,  is  taken  as  the  z  value.  The  z.  values  then 
range  Irom  0  to  I  on  each  crossplane.  The  advantage  of  the  mapping  method  is  that  equi  interval,  z  = 
constant,  coordinate  lines  arc  automatically  concentrated  in  the  region  of  large  curvature  where  the 
boundary  layer  characteristics  are  expected  to  vary  greatly. 

In  this  nonorthogonai  coordinate  system,  the  geodesic  curvatures  of  the  curves,  z  =  constant,  in 
the  direction  of  x  and  the  curves  \  -  constant  in  the  direction  of  z  are  given  in  Reference  7  as 
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and  the  geodesic  curvatures  of  the  curves,  y  =  constant,  in  the  directions  of  x  and  z  are 
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Since  y  is  selected  to  be  normal  to  the  hull  surface,  the  metric  coefficient  h3  can  be  chosen  as  unity. 
Along  the  x  and  z  directions,  hj  and  h-,  are,  in  general,  functions  of  x,y,  and  z;  that  is 


h,  =  h,  (x,  y,  z) 
h2  =  h2  (x,  y,  z) 


Expanding  h,  and  h2  in  series  about  the  coordinate  origin  leads  to 


hi  =  hio 


['”(i 


ah. 


3y  /y  =  0 


+  . . . 


(3 


h,  -  h 


20 


i  ah2' 


3y  h  =  0 


] 


where  h|0  and  h-,0  are  the  corresponding  values  at  the  surface  and  are  given  in  Reference  7, 


hio  _ 


+  +  /  3z 

\3x/  \3x/J 


1/2 


hi 


20 


'«Y  +  (ii 

\3z 


1/2 


(4 


By  comparing  Equation  (3)  to  Equation  (2)  and  dropping  the  higher  order  terms  in  the  series  expansion, 
we  can  write 
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a,  and  o,  are  defined  as 
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where  K130  and  K230  are  the  surface  curvatures  of  K13  and  K23.  Substitution  of  Equation  (5)  into 
Equation  (2)  results  in 
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The  other  related  geometric  parameters  are  given  in  Reference  7. 
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GOVERNING  EQUATIONS 


BOUNDARY  LAYER  EQUATIONS 

The  governing  second-order  boundary-layer  equations  for  a  steady,  incompressible,  laminar  or 
turbulent  flow  in  the  present  curvilinear  nonorthogonal  coordinate  system  are  given710  by: 
Continuity  Equation 
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x  Momentum  Equation 
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y-Momentum  Equation 


where  the  parameters  K,,  K2,  tx,  and  tz  are  defined  by 
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and  u,  v,  and  w  represent  the  velocity  components  in  the  x,  y,  and  z  directions,  respectively.  The 
parameters  p,  q  and  v  are  the  fluid  static  pressure,  density,  and  kinematic  viscosity;  gu'  v'  and  qv'  w'  are 
Reynolds  stresses.  The  non  slip  boundary  condition  is  applied  at  the  wall,  i.e.,  y  =  0.  At  the  edge  of  the 
boundary  layer,  i.e.,  y  =  <5,  Equations  (9)  to  (1 1)  can  be  written  as 
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Thus  at  y  =  3,  the  pressure  gradient  terms  can  be  expressed  in  terms  of  ue,  v£  and  we.  In  the  numerical 
examples  shown  later,  the  boundary-layer  Equations  (9),  (10),  (13),  and  (14)  are  solved  simultaneously 
through  the  procedures  described  in  the  following  sections. 

TURBULENCE  CLOSURE 
Cebeci  Smith  Eddy  Viscosity  Model 

Solutions  of  Equations  (9)  and  (10)  require  some  parameterization  of  the  shear  stresses  tx  and  tz. 
Cebeci7  used  an  isotropic  eddy  viscosity  concept  for  his  3-D  computation.  These  3-D  formulas  are  given 
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The  eddy  viscosity  tm  is  divided  into  the  inner  and  outer  layers.  They  are  continuously  connected 
at  y  yc.  where  (tm)j  =  (cm)0.  In  the  inner  region,  (cm)j  is  defined  as 
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where  K  is  a  mixing  length  approximation  equal  to  0.4  y  1 1  -exp(  1  -y/A)],  A  is  a  damping  factor  equal  to 
26v/ur.  and  uT  is  the  shear  velocity.  Also  y(r  is  an  intermittency  factor  to  account  for  the  turbulent  and 
nonlurbulent  flows  produced  by  the  fluctuating  outer  edge  of  the  boundary  layer,  and  is  defined  in 
Reference  11.  In  the  outer  region,  (£ni)0  becomes 

tem>»  =  00168  Y„  <1 


where  3(1*  is  the  planar  displacement  thickness  and  uIe  is  the  total  velocity  at  y  =  d. 

Modified  Eddy  Viscosity  Model 

In  order  to  calculate  a  thick  3  D  boundary  layer  by  means  of  the  eddy-viscosity  concept,  it  is 
necc'sary  to  include  the  curvature  effect1-  in  the  eddy  viscosity  formulations.  In  the  appendix,  a  formula 
is  dcmed  to  modify  the  inner  eddy  viscosity.  Equation  (17).  by 
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For  the  outer  layer  eddy  viscosity,  a  constant  value  0.0168  was  used  in  Cebeci’s  original  model,  as 
shown  in  Equation  (18).  This  has  been  proved  to  be  inadequate  for  a  thick  turbulent  boundary  layer13-14. 
A  suggested  revised  form  was  proposed  for  an  axisymmetric  case  in  Reference  14.  For  a  simple  extension 
of  Reference  14,  the  following  formulas  are  used  to  account  for  the  curvature  effect  on  the  outer-layer 
eddy  viscosity  when  K230  6  ^  0.2. 


<*Jo  =  F2  (0-0168)utedp  y(r 


where 


C  p  — 


^(1  +0.6  K230  d)2-I 
2.522  K22q  6 


\  he  value  of  C  F  becomes  \/(ro  +  0.6<5)2-ro2  12.5226  for  the  axisymmetric  thick  boundary  layer,  a  formula 
originally  proposed  by  Huang  et  al. 13,14 

Huang's  Eddy -Viscosity  Model 

From  three  sets  of  measurements1315,16  for  one  axisymmetric  and  two  3-D  bodies,  Huang  has 
obtained  a  single  correlation17  to  represent  the  eddy  viscosity  distribution  in  3-D  thick-stern  boundary 
lasers.  Huang's  turbulence  model  is  used  as  a  two-layer  isotropic  eddy-viscosity  model.  It  preserves  the 
(  ebeci  Smith  inner  layer  eddy  viscosity,  as  shown  in  Equation  (17).  The  outer-layer  eddy  viscosity  uses 
Equation  (18)  unless  Ap/A  ^  1.513  and  6Kj^q>  0.23,  where  Ap  is  the  effective  turbulence  area  defined 
in  Reference  17  as  60%  of  the  flow  passage  area  and  A  is  the  cross-sectional  area  at  the  same  axial 
location,  in  order  to  accommodate  the  present  numerical  method  described  in  the  later  section,  a 
constant  eddy  viscosity  selected  for  the  entire  outer  layer  is  given  by 
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where 


(£*'o  =  Max^  iTy)  +  (j;)  +2cosa  (j;)  (JiU  } 


i=  0.169  '/A^Tn  (y/c5)  exp[-  1.2(y/d)  -  1 .0667(y/c5)3] 

TRANSFORMED  BOUNDARY  LAYER  EQUATIONS 

The  boundary  layer  Equations  (8)  to  (10)  are  solved  in  the  form  of  transformed  equations  by  using 
the  following  transformed  variables: 


A  two  component  vector  potential,  which  satisifies  continuity.  Equation  (8),  is  defined  as 


y  =  /vs,ue  h20  sin  a  f  (4.  4,  rj) 
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where  u()  is  the  ship  speed.  Accordingly,  the  three  velocity  components  in  the  momentum  equations  can 
be  cast  as 
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Here,  primes  denote  differentiation  with  respect  to  rj.  and  all  the  coefficients  m,(i  =  1,2 . 27)  are 

defined  in  the  notation  section. 


In  order  to  further  reduce  the  complexity  of  the  boundary-layer  equations,  we  introduce  the 
following  parameters  to  replace  the  dependent  variables. 
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By  using  the  transformed  Equations  (22)  to  (25),  the  boundary-layer  equations  are  finally  written  as  five 
first  order  partial  differential  equations,  given  as: 
x  Momentum  Equation 
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Continuity  Equation 
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where  b  is  the  dimensionless  eddy  viscosity  and  is  equal  to  1  +£m/u.  The  nonslip  conditions  at  the  surface 
are  transformed  to  cf  =  w  =  G  =  0.  The  conditions  at  the  edge  of  the  boundary  layer  require  tT  =  1 
and  w  =  w/uQ.  Equations  (13)  and  (14)  for  the  pressure  gradients  at  the  boundary-layer  edge  are 
embedded  in  the  parameters  m,|  and  mp. 

EQUATIONS  ALONG  THE  KEEL  LINE 

In  Figure  1,  the  keel  line  is  identified  as  the  line  of  symmetry.  For  a  double  model,  the  water  line 
also  becomes  a  symmetry  line.  On  this  symmetry  plane,  w  and  0p/0z  are  zero,  thus  causing  the 
z  momentum  equation  to  be  singular.  However,  differentiation  with  respect  to  z  yields  a  nonsingular 
equation.  After  performing  the  necessary  differentiation  for  the  z-momentum  equation  and  applying  the 
appropriate  symmetry  conditions,  we  obtain  the  following  equations: 

Continuity  Equation 
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This  set  of  equations  can  be  transformed  by  using  a  similar  set  of  transformed  parameters 
described  in  Reference  7,  with  the  same  modifications  as  used  in  the  previous  section.  The  final 
transformed  first-order  partial  differential  equations  can  be  summarized  as  follows: 
x  Momentum  Equation 
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Continuity  Equation 
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The  condition  at  the  edge  of  the  boundary  layer  is  changed  tow  =  w/c/u().  All  the  other 
parameters  are  the  same  as  the  ones  for  the  general  3-D  case  except: 
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NUMERICAL  METHOD 
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NUMERICAL  CALCULATION 


Since  the  governing  boundary-layer  equations  derived  previously  are  parabolic,  a  numerical 
marching  scheme  is  adequate  for  obtaining  a  solution.  The  marching  process  is  started  at  an  initial 
station  in  the  bow  region  and  is  continued  downstream.  In  the  crossflow  direction,  the  calculation 
proceeds  from  the  ked  line  to  the  water  line.  The  flow  at  the  initial  station  was  solved  by  assuming  all 
m,  =  0  except  m(  =  0.5.  which  is  a  flat  plate  solution  Equations  (34)  to  (36)  along  the  keel  line  are 
then  employed  to  start  the  cross-stream  sweep.  Kclle'  .■>  box  method7  is  used  as  the  numerical  solution 
scheme.  This  method  is  accomplished  by  first  differencing  th^  transformed  boundary-layer  Equations  (26) 
to  (30)  with  central  differences.  Newton’s  method  is  then  applied  to  linearize  the  difference  equations 
and,  finally,  the  linearized  system  is  solved  bv  a  block  elimination  method.  Details  of  this  scheme  may  be 
found  in  Reference  18. 

The  computations  were  carried  out  -.n  the  DTNSRIX  CYBER  176.  Typical  run  times  were  80  s 
for  one  boundary-layer  sweep  on  an  87  x  1 5  surface  grid  with  an  average  of  60  normal  grid  points. 


PR(X  EDURE  FOR  UPDATING  OUTER  EDGE  BOUNDARY  CONDITIONS 

An  iterative  procedure  between  the  boundary  layer  solution  and  the  potential  solution  is  generally 
necessary  in  order  to  account  for  the  strong  interaction  between  these  two  regions  of  ship  stern 
flows.--6’8  This  interaction  is  usually  considered  by  using  the  displacement  body  concept  at  y  d*  or  the 
equivalent-source  concept  at  y  -  0  (Reference  I  hi.  The  disadvantage  of  using  the  former  approach  is 


that  a  new  grid  of  the  displaced  body  must  be  generated  for  each  iteration.  Although  the  latter  approach 
avoids  the  complexity  of  regenerating  the  grid  for  each  iteration,  a  divergent  source  distribution  is  often 
encountered  at  the  stern  region.  Improper  treatment  of  this  ill-behaved  source  distribution  causes  the 
global  iteration  to  be  divergent.20,21  This  problem  becomes  more  serious  when  flow  separation  exists  in 
the  stern  region. 

In  order  to  validate  the  present  program  for  3-D  curvature  effects  and  to  avoid  the  problems 
associated  with  the  displacement-body  and  equivalent-source  approaches,  a  simple  procedure  was  used  to 
account  for  the  thick  boundary-layer  effect.  In  this  procedure,  firstly  (i)  the  boundary  conditions  at  y  =  6 
are  obtained  from  the  potential-flow  solution  at  the  body  surface  for  the  first  iteration.  For  the  second  or 
higher  iterations,  these  conditions  are  recalculated  as  the  off-body  points  from  the  singularity  distribution 
for  the  original  body.  The  locations  for  applying  the  boundary  conditions  are  obtained  from  the  previous 
boundary-layer  solutions.  Secondly,  (ii)  the  normal  pressure  gradient,  as  given  by  Equation  (111,  is  solved 
by  varying  dp/dx  and  dpldz  linearly  from  y  =  0  to  y  =  6  for  the  curvature  effect  only. 

COMPARISON  OF  EXPERIMENTAL  DATA  WITH  NUMERICAL  PREDICTIONS 

Wind-tunnel  test  data13,15,16  have  been  collected  at  DTNSRDC  for  the  so-called  axisymmetric 
Afterbody  1  and  for  two  three-dimensional  bodies  of  2:1  and  3:1  elliptic  cross-sections.  The  latter  two 
three-dimensional  bodies  are  generated  from  axisymmetric  Afterbody  1  by 

y2  -2 

—  +  3  z‘  =  r0M  (37) 

u 

for  an  elliptic  cross-section  on  the  yz-plane;  r0  is  the  corresponding  axisymmetric  radius  at  x,  and  a  is 
equal  to  2  or  3,  respectively,  for  the  2:1  or  3:1  body.  A  plot  of  cross-sections  at  the  same  x  location  for 
these  three  bodies  is  shown  in  Figure  2.  All  three  models  have  a  length  L  =  3.066  m  and  were  tested  in 
a  wind  tunnel  at  a  nominal  speed  of  30.48  m/s.  Figure  3  shows  the  surface  longitudinal  and  transverse 
curvatures  of  Afterbody  1  in  the  stern  region.  Since  the  order  of  magnitude  of  the  longitudinal 
curvatures  for  the  two  3-D  bodies  is  the  same  as  that  of  Afterbody  1,  Figure  3  only  shows,  at  a  typical 
station,  the  transverse  curvatures  K2 -j0  'n  tlic  girthwise  direction,  indicated  by  9  as  defined  in  Figure  2, 
for  both  3  D  bodies.  The  3:1  body  has  the  least  transverse  curvature  when  6  is  between  0°  and  80°.  It 
also  possesses  the  largest  K -> 30  value  when  9  exceeds  80°.  and  reaches  a  maximum  at  9  -  90°. 

Boundary  layer  computations  were  done  first  for  the  flow  around  axisymmetric  Afterbody  I  in 
order  to  validate  the  ability  of  the  program  to  include  curvature  effects  across  the  boundary  layer.  The 
computations  were  started  with  laminar  flow  and  transitioned  to  turbulent  flow  at  x/L  =  0.03.  The 


pressure  coefficient  distribution  computed  from  the  potential  flow  solution  for  the  displaced  body14  is 
plotted  in  Figure  4a  as  compared  to  the  measured  distribution.  This  calculated  pressure  distribution  was 
used  as  the  boundary  condition  for  the  following  program  validation  runs.  In  Figures  4a  and  4b,  five  sets 
of  velocity  predictions  are  shown  along  with  the  measured  velocity  profiles.  These  calculations  are  (i!  3  D 
thin  boundary-layer  calculation,  indicated  by  K)3  =  K23  =  0  in  Figure  4a;  (ii)  3-D  boundary-layer 
calculation  with  only  transverse  curvature  effect,  indicated  by  K)3  =  0,  K->3  in  Figure  4a;  (iiil 
axisymmetric  boundary-layer  calculation2-  including  only  the  transverse  curvature,  indicated  by 
axisymmetric  calculation  in  Figure  4b;  (iv|  3-D  boundary-layer  calculation  with  both  longitudinal  and 
transverse  curvature  effects,  indicated  by  Kj3,  K.-,3  in  Figure  4a;  all  the  above  four  calculations  use 

Cebeci-Smith  turbulence  closure,  (v)  3-D  boundary-layer  calculation  with  both  curvature  effects  and  the 
modified  turbulence  model,  indicated  by  K|3,  K73,  and  modified  turbulence  modeling  in  Figure  4a.  The 
axisymmetric  boundary-layer  calculation  (iii)  solves  the  boundary-layer  equations  and  includes  the 
transverse  curvature  effect  for  axisymmetric  flow.  In  terms  of  numerics,  calculation  (iii)  is  an  exact 
solution  of  calculation  (ii).  Comparisons  of  calculation  methods  (ii)  and  (iii)  are  shown  in  Figure  4b 
Velocity  predictions  from  all  calculation  methods,  except  (iii),  are  shown  in  Figure  4a.  As  suggested  by  its 
order  of  magnitude,  the  transverse  curvature  exerts  a  significant  influence  on  the  velocity  predictions, 
particularly  in  the  outer  region  of  the  viscous  layer.  The  effect  of  the  longitudinal  curvature  on  the 
velocity  prediction  shows  no  influence  at  x/L  =  0.846  and  1  to  ~  2%  in  the  region  close  to  the  wall  at 
x/L  =  0.934.  More  importantly  in  this  “inner  region”,  the  velocity  prediction  in  Figure  4a  using  two 
different  turbulence  models  (Cebeci  Smith  and  modified)  are  significantly  different.  The  good  agreement 
of  calculations  (ii)  and  (iii)  indicates  the  need  to  include  the  curvature  effects  in  3  D  boundary-layer 
computations.  The  good  agreement  between  the  predictions  of  calculation  (v)  and  the  measured  velocity 
profiles  suggests  that  the  modified  turbulence  mode!  gives  reasonable  predictions  for  axisymmetric  flow. 

Based  on  the  findings  for  axisymmetric  Afterbody  I,  predictions  of  the  flow  about  the  two  3  D 
bodies  were  made  for  the  following  approximations  in  order  to  examine  the  significance  of  the  curvature 
effects  and  the  turbulence  modeling.  Thus,  they  are 

(i)  Thin  boundary-layer  calculation  with  the  Cebeci-Smith  model  and  using  the  calculated  potential 
flow  solution  at  the  wall  1C  I  as  the  boundary  conditions, 
liil  Thin  boundary  layer  calculation  with  the  Cebcci  Smith  model  and  using  the  boundary 
conditions  1C  )  obtained  from  the  potential  solution  at  the  edge  of  the  boundary  layer. 

The  following  three  calculations  use  C  and  the  calculated  values  of  K(1  and  K->3  obtained  from 
liquation  i7).  The  boundary  layer  calculations  were  performed  by  using  the  procedures  described 
previously  for  updating  the  edge  boundary  conditions.  The  turbulence  model  for  each  calculation  is: 

(iii)  Cebeci  Smith  turbulence  model 
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(iv)  Modified  turbulence  model 

(v)  Huang’s  turbulence  model  (v) 

The  potential  flow  solution  was  obtained  by  using  the  XYZ  potential  flow  code,23  which  utilizes  Hess 
Smith’s  method24.  Calculation  methods  (i)  and  (ii)  show  the  effect  of  the  thick  boundary  layer,  and 
calculation  methods  (ii)  and  (iii)  show  the  effect  of  the  added  curvature  terms.  The  results  of  the  last 
three  calculations  illustrate  the  significance  of  the  turbulence  modeling. 

Figure  5  shows  plots  of  two  different  calculated  Cp  distributions  used  as  the  boundary  conditions 
Cpw  and  C  .  Measured  values  of  Cpw  and  are  also  shown  in  Figure  5  for  comparison.  Although  the 
measured  data  vary  from  station  to  station,  this  figure  indicates  that  the  pressure  variation  across  the 
boundary  layer  becomes  large  when  0  increases.  The  calculated  C__  values  agree  reasonably  well  with  the 

pc 

measured  C  values.  The  calculated  Cn„.  values,  however,  have  too  strong  a  girth-wise  adverse  pressure 

pe  pw 

gradient  for  0  between  80°  and  90°.  This  behavior  is  responsible  for  the  earlier  predicted  separation, 
which  will  be  discussed  later.  The  predicted  drop-off  of  values  near  6  =  90°  and  x/L  >  0.89  is  also 
associated  with  this  earlier  predicted  separation.  With  the  CnUiand  C..  values  shown  in  Figure  5,  the 

pw  pc 

actual  magnitudes  of  the  crossflow  velocities  used  in  the  computations  introduced  in  Equations  (13)  and 
(14)  are  plotted  in  Figure  6.  Although  the  experimental  data  are  quite  scattered,  the  predicted  we  values 
are  somewhat  larger  than  the  measured  values,  except  at  x/L  =  0.719  for  the  3:1  body. 

The  predicted  axial  velocity  profiles  are  shown  in  Figures  7a  and  7b  for  the  2:1  body  and  Figures 
8a  and  8b  for  the  3:1  body.  Results  of  calculation  (ii)  are  shown  in  all  figures  as  baselines  for  comparing 
the  relative  magnitudes  of  the  predictions.  Separation  is  predicted  when  flow  reversal  occurs  in  the  axial 
direction  and  the  calculation  stops  at  the  predicted  separation  point.  Therefore,  at  each  0-value,  no  profile 
is  plotted  downstream  of  the  separation  point.  Table  1  shows  the  predicted  separation  locations  x  and 
0sep  for  each  case  as  compared  to  the  experimentally  observed  locations  of  separation  for  the  2: 1  body 
and  Table  2  for  the  3:1  body.  Calculation  set  (i)  predicts  early  separation.  The  other  four  calculation  sets 
predict  the  locations  of  separation  inception  within  3%  of  the  observed  values.  Figure  3  indicates  that  at 
x/L  =  0.894  the  transverse  curvature  for  axisymmetric  Afterbody  I  equals  those  of  the  two  3  D  bodies 
near  0  =  80°.  As  0  increases  further,  the  curvature  becomes  larger  for  the  two  3  D  bodies.  But  the 
predicted  velocity  profiles  in  Figure  4  for  axisymmetric  Afterbody  1  respond  more  strongly  to  the 
curvature  effect  than  those  in  Figures  7  and  8  for  the  two  3  D  bodies.  This  is  due  to  the  three 
dimensionality  of  the  latter  predictions. 


Figure  8  (Continued) 
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TABLE  1  -  PREDICTED  LOCATIONS  OF  SEPARATION  INCEPTION  COMPARED  TO 


THE  EXPERIMENTAL  OBSERVATION  FOR  THE  2:1  BODY 


Boundary 

Curvature 

Turbulence 

<*/L>sep 

0sep 

Condition 

Effects 

Modeling 

Experiment 

0.90 

80°  ~  90° 

Predictions 

(i)  Wall 

No 

CS* 

0.775 

90° 

lii)  Edge 

No 

CS 

0  894 

80.7° 

(iii)  Edge 

Yes 

cs 

0.934 

83.2° 

(iv)  Edge 

Yes 

Modified 

0.914 

85.6° 

(vi  Edge 

Yes 

Huang 

0914 

8C).9° 

*Cebeci  Smith 


TABLE  2  -  PREDICTED  LOCATIONS  OF  SEPARATION  INCEPTION  COMPARED  TO  THE 
EXPERIMENTAL  OBSERVATION  FOR  THE  3:1  BODY 


(x/l'sep 


Boundary 

Curvature 

Turbulence 

Condition 

Effects 

Modeling 

Experiment 

Predictions 

(ii  Wall 

No 

CS* 

litl  Edge 

No 

CS 

(nil  Edge 

Yes 

CS 

it  vi  Edge 

Yes 

Modified 

(vi  Edge 

Yes 

Huang 

(.834 

90° 

1.934 

75° 

1.956 

75° 

1.934 

85.6° 

1.907 

90° 

Cehcei  Smith 


The  last  three  calculation  sets  for  both  3  1)  bodies  show  the  effect  of  turbulence  modeling.  Huang's 
turbulence  model  (calculation  method  ivii  underestimates  the  inner  region  velocity  as  compared  to  the 
Cebeci  Smith  model  (calculation  method  Him.  Because  of  the  overprediction  of?,  in  Equation  (21l  for  the 
modified  turbulence  model,  calculation  method  (tvi  gives  the  worse  predictions.  This  overprediction  of  t  ( 
is  related  to  the  use  of  local  curvature  and  boundarv  (aver  parameters.  Although  the  Ceheci  Smith 
turbulence  model  does  not  include  curvature  effects,  it  works  better  in  terms  of  velocity  prediction  than 
the  other  two  turbulence  models  for  the  two  tested  3  I)  bodies.  Figure  9  shows  the  eddv  viseosilv 
comparisons  fetween  the  measured  data  and  the  calculated  values  by  the  Cebeci  Smith  model  and 


Huang's  turbulence  model  (indicated  by  dashed  line).  Huang's  turbulence  model  closely  predicts  the  levels 
of  the  eddy  viscosity  in  thick  boundary  layer,  as  compared  to  the  measured  data,  but  it  underpredicts  the 
velocity  profile  in  the  inner  region  unlike  the  Cebeci  Smith  model. 

Generally,  the  predictions  of  calculation  method  (iiil.  including  the  ones  at  x/L  =  0.719  and  0.810. 
and  6  =  90°  for  the  3:1  body,  are  in  good  agreement  with  the  measurements.  The  deviation  of  the 
predicted  velocity  from  the  measured  velocity  at  0  =  90°  and  x/L.  =  0.854:  0.894  for  the  3:1  body  in 
Figure  8  may  be  related  to  the  proximity  of  the  measurements  to  the  separation  point. 

CONCLUSION 

Three-dimensional  thick  boundary  layer  equations,  including  variations  of  the  metric  coefficients, 
surface  curvatures  and  static  pressure  coefficients  with  the  normal  distance  from  the  wall,  are  solved  by 
using  the  marching  technique  of  Keller's  box  method.  A  nonorthogonal  surface  fitted  coordinate  system  is 
used.  Formulas  for  calculating  second  order  geometric  parameters  are  derived.  Three  turbulence  models, 
two  of  which  include  the  curvature  effects,  using  the  eddy  viscosity  concept  are  compared.  Three  sets  of 
measured  data,  for  one  axisymmetric  and  two  3-D  flows  are  compared  with  calculated  results.  The 
conclusions  drawn  from  the  present  study  can  be  summarized  as  follows: 

(i)  The  3-D  boundary  layer  computations  with  varying  geometric  curvatures  in  the  normal 
direction  have  been  validated  by  comparing  the  axisymmetric  boundary  layer  calculations  to 
experimental  data. 

Iiil  The  predictions  show  a  significant  improvement  in  the  present  solution  procedures  by  including 
the  variation  of  transverse  and  longitudinal  curvatures  across  the  thick  turbulent  stern  flows. 

liii)  The  calculation  procedure  that  uses  the  edge  boundary  conditions  (C  t  for  a  thick  boundary 
layer  predicts  a  more  accurate  flow  separation  location  for  a  separated  How  than  that  using  the 
surface  boundary  condition  !('  I. 

liv)  The  modified  turbulence  modeling,  which  uses  the  local  curvature  and  local  boundary  layer 
parameters,  works  reasonably  well  for  axisymmetric  flows  but  not  as  well  for  3  D  flow. 

Ivi  Although  Huang's  turbulence  model  correlates  eddy  viscosity  distributions  in  thick  boundary 
layers  better  than  the  C  ebeci  Smith  model,  the  velocity  predictions  using  Huang's  model 
underestimate  the  inner  region  velocity . 
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APPENDIX 

DERIVATION  OF  INNER  LAYER  EDDY  VISCOSITY  E Q.  (19) 

For  a  small  pressure  gradient  flow,  the  x-  and  z-monientum  Equations  (9)  and  (10)  in  the  inner 
region  of  the  viscous  layer  can  be  approximated  by 


d  d 

—  (h,  h-,  r.l  =  —  (h ,  h-,  t  l  =  0 

3y  1  -  x  3y 


Integrating  Equation  (A.I)  leads  to 

ai  a2  Ts  =  Ts«  =  e  u?  <A 

The  subscripts  s  and  w  stand  for  the  streamline  direction  and  the  wall.  By  introducing  the  following 
nondimensional  quantities. 


u*  .  v 


^  =  foL 
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v  v 


^  no  ~ 


K  nn  v 


k23()  " 


K-n,,  v 


°\  =  °l 


<*->  =  a*, 


and  by  means  of  the  eddy  viscosity  concept,  the  Equation  (A. 2).  in  the  fully  turbulent  flow,  can  be 
written  in  the  nondimensional  form  as 


du  *  du  ‘ 

_  ^  *  _ 

dY  *  t'"  dY7 


n  here 


.  In  la -jo .*  i 

/  - =_J —  .  a;  *  0  * 


By  comparing  Equation  (A.4)  to  the  formulation  of  the  eddy  viscosity  as  shown  in  Reference  12  for  2  D 
flows.  Equation  (19)  can  be  obtained. 
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